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We study electronic transport in junctions consisting of a superconductor electrode and two ferro- 
magnet (F) leads in which crossed Andreev reflections (CAR) and elastic co-tunnelings are accom- 
modated. We model the system using an extended Blonder-Tinkham-Klapwijk (BTK) treatment 
with a key modification that accounts for spin-dependent interfacial barriers (SDIB). We compute 
current-voltage relations as a function of parameters characterizing the SDIB, magnetization in the 
F leads, geometry of the junction and temperature. Our results reveal a rich range of significantly 
altered physics due to a combination of interfering spin-dependent scattering processes and popu- 
lation imbalance in the ferromagnets, such as a significant enhancement in CAR current and a sign 
change in the relative resistance between the parallel and antiparallel alignments of the magnetiza- 
tion of the F leads. Our model accounts for the surprising experimental findings of positive relative 
resistance by Colci et al. [Phys. Rev. B 85, 180512(R) (2012)] as well as previously measured 
negative relative resistance results, both within sufficiently large parameter regions. 

PACS numbers: 74.45. +c, 74.25. F-, 74.25.fc, 72.25.-b 



I. INTRODUCTION 

Over the past decades, extensive theoretical and ex- 
perimental research has focused on electronic transport 
in superconductor-normal metal (S-N) heterostructures, 
wherein conducting electrons propagate through bulk 
materials and scatter at interfacea^Hl. More recently, at- 
tention has turned to systems comprising an S electrode 
in contact with multi-N terminals, in which a unique scat- 
tering process known as the crossed Andreev reflection 
(CAR) can occui^EI The CAR process involves an elec- 
tron in an N terminal impinging the S accompanied by 
a hole of opposite spin reflecting in another N terminal 
separated within the superconducting coherence length, 
and generating a Cooper pair in S . The non-local and 
coherent nature of CAR makes such a device a platform 
for exploring quantum entanglement, providing potential 
application to quantum computinglS 

Another avenue of S-N junction studies which has re- 
cently gained prominence is the presence of ferromagnetic 
order in the normal system. Transport physics can be 
significantly altered when the normal system is a ferro- 
magnet (F). For instance, the F bulk can be created by 
spin imbalance, changing the populations involved in the 
transport processes as well as spin-dependent interfaces 
that modify both the proximity effedP*"IH anc | electron 
scattering-ISJES j n g_p junctions, the interplay between 
ferromagnetic and superconducting orders around the in- 
terface strongly affects Andreev reflections^^!] (the di- 
rect process of an incoming electron reflecting as a hole of 
opposite spin within the same terminal 2 ^) , while in S-F-S 
systems it can lead to Josephson 7r-j unctions sustaining 
negative critical currents^HSH. Recent experiments^^ 
have been performed on S-FF-S devices in which two F 
bridges are both laid across two S electrodes, with the 
separation between the two F bridges being smaller than 




FIG. 1: (Color online) Cartoon representation of an SFF junc- 
tion. The system has two ferro-magnetic leads Fl (F2), illus- 
trated as the shaded regions in the second (third) quadrant 
having width Wjp) and separation L between them, and a su- 
perconducting electrode, S, illustrated as the shaded region 
in the x > plane with width Ws- The Fl(2)-S interfaces 
(twilled regions) are described by spin-dependent interface pa- 
rameters Z\ and Zi, respectively. Black dashed arrows indi- 
cate six scattering processes experienced by an electron in Fl 
incident on the interface (Inc., white dashed arrow): local nor- 
mal reflection (LNR) , local Andreev reflection (LAR) , crossed 
normal reflection (CNR), crossed Andreev reflection (CAR), 
quasi-particle transmission (QPT) and quasi-holc transmis- 
sion (QHT). 



the superconducting coherence length (thus capable of 
accommodating CAR) . The re sults show that the system 
can exhibit both positive^^ and negative^ relative re- 
sistances, SR, between two different cases characterized 
by an anti-parallel (AP) or parallel (P) alignment of mag- 
netization in the F leads, respectively. Since coherent 
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transport is not revealed in the experimental data, the 
S-FF-S device can be considered as a series connection of 
two independent S-FF junctions (as illustrated in Fig.[T]). 
In such cases, the negative SR can be understood by 
an intuitive picture that the suppression of CAR due to 
spin imbalance raises the resistance of the P case, or can 
be explicitly explained by a Blonder-Tinkham-Klapwijk 
(BTK) treatment with spin-independent interfacial bar- 
riers^. However, such a BTK model can not explain the 
counter-intuitive data of positive which suggests that 
a competitive effect should be incorporated. 

Motivated by these experiments, we provide one of the 
first comprehensive theoretical studies of S-FF junctions 
having spin-dependent interfacial barriers (SDIB). The 
presence of SDIB at the junction between a supercon- 
ductor and a single ferromagnet has been shown to con- 
siderably alter transport properties^* 2 !!. For our situ- 
ation, the SDIB plays a prominent role in affecting co- 
herence and cross correlations between excitations in the 
two ferromagnets. We explore the large parameter space 
of the system and show that SDIB can give rise to a 
rich range of physics depending on the choice of param- 
eters. We provide a consistent scenario explaining both 
competing experimental results described above. In fact, 
a small subsection of the studies investigated here has 
been directly employed to explain the counterintuitive 
experimental results presented in Ref. [35] We propose 
a microscopic picture describing the SDIB within an ex- 
tended BTK model and compute transport properties of 
the system as a function of the derived SDIB parameters. 
We carefully discuss the various scattering processes re- 
sponsible for the difference in behavior between the P and 
AP cases and show that interference between these pro- 
cesses can play a significant role. We analyze transport 
properties as a function of magnetization, the geometry 
of the system, and temperature. We show that they re- 
veal several interesting features, such as, the sought-after 
differing behavior for the P and AP cases, oscillations as 
a function of geometry, and, under certain conditions, 
enhancement of crossed Andreev reflection, which would 
be indicative of long-range entanglement. 

Our studies thus not only reconcile the apparently con- 
flicting results of the various experiments, they also pro- 
vide several other aspects to explore within the same ex- 
perimental settings. Our extended BTK model also ap- 
plies to several other possible experimental geometries, 
including other multi-terminal superconducting hybrids 
that possess a two dimensional nature characterizing the 
elongation and separat ion o f the terminals, such as multi- 
layeitEEIH or multi-wir d 39 * 40 ! devices in the presence of the 
spin-dependent interfacial scattering. 

The paper is outlined as follows. In Sec. [n] we model 
the geometry of the system and present the correspond- 
ing Hamiltonian describing six scattering processes. We 
also discuss a possible microscopic mechanism resulting 
in an SDIB. In Sec. |III[ by applying the BTK treatment 
we derive the current-voltage (/-V) relation as a func- 
tion of exchange energy, SDIB parameters, the geometry 



of the junction and temperature. In Sec. IV we illus- 



trate data showing the dependence of transport proper- 
ties on several groups of variables associated with fer- 
romagnetism, geometry of the system and temperature, 
respectively. We discuss the effects of SDIB compared 
to the spin-independent case and other effects competing 
or cooperating with SDIB. In Sec. |V]we summarize our 
results and comment on the scope of our work. 



II. THE MODEL 

To model our basic set-up shown in Fig. [T] we consider 
a two dimensional junction consisting of an S electrode 
in contact with two F leads. The S electrode of width 
Ws is located in the x > half plane, while the two F 
leads (Fl and F2), of widths W\ and W2 respectively, are 
located in the x < half plane and separated by L in the 
y direction. This system is described by a Bogoliubov-de 
Gennes (BdG) Hamiltonian 



H - ae ex (r) + ^(r) 



A(r) 



A(r) 

[H - ae ox (r) + l/i(r)] 



■(1) 



Here Ho = — — €f is the free Hamiltonian with 

Fermi energy, ep, r = (x,y) denotes the coordinates, 
a = ± denotes up (down) spin states, respectively, and 
a = — a. The superconducting gap A, exchange energy 
e cx and interface potential Vj are given by 



A(r) = A0(x)d s (y) 
6 ox (r) = 6{-x) J2 <?Wv) 



(2) 
(3) 



3=1.2 



Vj(r) = S(x) J2 Z MV), ( 4 ) 

3=1,2 

where e| x j 2 are exchange energies in Fl(2), respectively, 
indicating the Zeeman splitting between up and down 
spins induced by ferromagnetism. The parameter Zj a 
describes the corresponding S-F interfacial barrier, 9 is 
the step function, and 
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Here the 5 function and the step function are applied 
to confine the physical quantities in the corresponding 
regions. 

The SDIB parameters Z a can be decomposed as 



Zrr — Zn i Z s 



(7) 



where the spin-independent component, Zo, represents 
effects of an oxide layer or the local disorder at the inter- 
face^, while the spin-dependent one, Z s , can be induced 
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by ferromagnetism. The sign depends on whether a par- 
ticular spin component a is the majority (+) or minority 
(— ) carrier. As one of the prevalent causes for a SDIB, we 
propose a semiclassical picture to derive Z s as a function 
of the average variance of exchange energy, (Se cx ) , over 
a microscopic length, £, associated with the interfacial 
properties: 



= tan 



(fc ex )£ 



(8) 



This form can be derived by considering the phase accu- 
mulation of an electron in a microscopic model involving 
propagation through a finite layer having spin-dependent 
potentials and relating it to the phase shift in an effec- 
tive model described by the BTK ((-function barrier^ 
(see Appendix [AJ. Equations ^ and ^ show that for 
a purely magnetic interface (Zq = 0) one has Z^ = —Z^ 
with the magnitude sensitive to the microscopic configu- 
ration. 

Now we turn to six different scattering processes 
through the interface, which describe the transport 
physics in the system. As illustrated in Fig. [IJ given an 
electron of energy lower than the superconducting gap in- 
jected from Fl, there are two local scattering processes at 
the Fl-S interface, (1) the direct backscattering, or local 
normal reflection (LNR), and (2) local Andreev reflection 
(LAR), as well as two crossed scattering processes medi- 
ated by the superconducting order at the F2-S interface, 
(3) electron backscattering called elastic co-tunncling, 
and (4) crossed Andreev reflection (CAR). For conve- 
nience in comparison, we refer to the elastic co-tunneling 
as crossed normal reflection (CNR) in what follows. If 
the energy is higher than the superconducting gap, we 
have two more processes: (5) quasi-particle transmissions 
(QPT) and (6) quasi-hole transmissions (QHT) into the 
S region. In Sec. Ill we apply the BTK treatment to 
compute scattering amplitudes subject to the incoming 
state as well as the SDIB and hence obtain a net charge 
current carried by these scattering processes. 



III. BTK TREATMENT 

The BTK treatment is to consider charge transport as 
a net effect of reflections and transmissions of electrons 
or holes at the N-S (F-S) interface. The scattering am- 
plitudes of reflections and transmissions are obtained by 
solving the BdG equation with the interface potential. 
In the presence of voltage drop across the interface, an 
induced current is computed by summing the probability 
current contributed by each scattering process, weighted 
by the Fermi distribution. The calculations are aimed 
at a current-voltage relation as a function of physical 
variables of interests, such as the interface parameters, 
exchange energy, geometry of the system, and tempera- 
ture. 

We start with the Hamiltonian of Eq. for con- 
sidering incoming and outgoing waves incorporating the 



six scattering processes (discussed in Sec. As pre- 
scribed by the BTK treatment, we solve for the forms 
of the incoming and outgoing waves by boundary con- 
dition matching. Our theoretical set-up requires careful 
accounting of the spin-species, the multiple channels and 
the SDIB; we thus provide a detailed outline of the pro- 
cedure below. 

We consider an incoming wave of energy E from the 
F side (x < 0). The wave function is of the form of 
a plane wave in the x direction multiplied by a bound 
wave in the y direction. We use indexes t, a indicating 
a particle (r = +) with spin a or a hole (r = — ) with 
spin (7, the channel number n labeling the bound states 
in the y direction, and an index j denoting the wave in 
Fj (j = 1,2) regions (see Fig. [T]). The incoming wave is 
written as 



rajn 



5 T+ 






Sr- 



(9) 



where S is the <5-function. The wave vectors, p, and the y 
component wave function, $, are given below in Eqs. ( 12 ) 



and (14), respectively. The sign in front of of the wave 



vector is chosen to match the direction of the group ve- 
locity. 

The outgoing wave is represented as a linear combina- 
tion of degenerate scattering modes with energies and 
group velocities corresponding to the incoming wave. 
The wave functions \&j/s(r) in Fj/S regions are sepa- 
rately given as 
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i=i 

Mi 
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*?(y) 



i=i 



E 6 H I )e- i ^Ui(y), 



E 

i=i l 



(10) 



sf(y). 
(ii) 



Here a\ is the amplitudes for LAR (CAR) of channel I if 
j is the same as (different from) the incoming wave, and 
similarly bf represents the LNR (CNR) processes. The 
amplitudes q (di) correspond to quasi-particle (quasi- 
hole) transmissions in channel I. These amplitudes of the 
outgoing waves are also functions of the indices {t, a, j, n} 
of the incoming wave, which have been dropped here for 
convenience. The quasi-particle basis u and v satisfy 
u 2 = 1 - v 2 = |(1 + y/l - A 2 /E 2 ). The wave vectors of 
the x component of the wave functions in channel I are 
given by 



Pl'l = Jl + rE + aef-iln/W,) 



(12) 



1 + t\J E 2 - A 2 - (ln/W s ) 2 . (13) 
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From here on, we take the Fermi energy and inverse 
of the Fermi wave vector as energy a nd length units 
(e F = fcp 1 = 1), respectively. In Eq. (10), the upper 



bounds of the summations, M^ b , is given by the highest 
current-carrying mode, above which a mode has a purely 
imaginary wave vector and hence carries no current. In 
Eq. ( 11 ), the wave vector can never be purely imaginary 
if E < A. In such case, we choose Ms large enough to 
guarantee the convergence in numerical calculations^. 
The y components of the wave functions in Eq. ( 12 ) and 
( 13 1 are given by 



and project it onto the channel m in the Fi region, 
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(15) 



To solve for the amplitudes, we match the incoming 
and outgoing waves by imposing the boundary conditions 
at the interface (x = 0): 



Using Eq. ( 18 ) to replace q and di in Eq. ( 20 1 , we rewrite 



the right-hand side of Eq. ( 20 1 as 
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where fltl„ is a 2 x 2 matrix defined as 



,(21) 



in (o, y) + J2 *° ut (o, y) = n ut (o, y), 
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(17) 



Here Eq. (|16[) is the continuity equation, while Eq. (17 1 is 



obtained from integrating the BdG equati on through the 

interface (J °_ dx). For Zi„ = Z%&, Eq. jl7j ) reduces to 
the boundary condition for the case of a spin-independent 
interface discussed in Ref. 1351 We substitute the wave 
functions in Eqs. (|9 11 ) into Eq. ( 16 ) and project it onto 
channel m in the S region, 



3t+ | v 



E 



E 

i=i 



m i 

I )E A ^ 
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(18) 



where 



A J - 



dy#f'(y)^(l/). 



(19) 



Similarly, we substitute the wave functions into Eq. (17) 



ST- 7 = 



U V 
V u 



M S 



i=i 



(22) 



Combining Eqs. (20 22) and letting i run from 1 to 2 as 



well as m run through all current-carrying channels, we 
get a set of linear equations for obtaining the amplitudes 
a\n and b l m . 

Next we turn to compute the probability currents, 
which in turn yield the current-voltage relation. The 
probability currents for wave functions of a Nambu form, 

( ipi ip2 ) j are defined as 



Im(^*Wi - $jV^2), 



(23) 



where the ipi term represents the contribution from elec- 
trons and the ip2 term does that of holes. The wave func- 
tions in Eqs. (plTTO) carry no current in the y direction. 



The outgoing probability currents in the x direction at 
the interface (x = 0) in the ferromagnet and supercon- 
ductor regions [J(E) and J (E), respectively] are given 
by 



rajn 



h 
m 



5 T ,_ 



Re(?V'i') 
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x(n 2 -H 2 )/./ : 



rajni 



(25) 
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normalized by the incoming current, 

h 



J" 



(26) 



Here the subscripts (superscripts) denote the correspond- 
ing incoming (outgoing) state. In Eq. (24 1, {r',cr'} = 



±{t, a} represents normal or Andreev reflections, respec- 
tively, while j' = j or j (the counter part of j) denotes 
local or crossed processes, respectively. 

We assume the bias voltage, V, to be the same across 
both F leads and the S electrode, as set up in the exper- 
iment of Ref. [35] We remark that our method can also 
be applied to the case of different voltages on different 
leads. Such a scenario could enable one to probe other 
physical properties, such as current correlation and long 
range entanglement, in the system. 

Following the standard BTK formalisrrP^U, we obtain 
the charge current, /, carried by incoming particles or 
holes (r = ±) with a spin in channel n in lead Fj (Ap- 
pendix [B]) , 

oo 

I TtTjn (V) = re J dE 
o 

x {C 1 - E E\ J %n'\)[fo(E - reV) - f (E)] 

j' n' 

+ EE \j^n'\ifo(E) f (E + reV)}}, (27) 

j' n' 

up to a constant associated with density of states, Fermi 
velocity and an effective cross-sectional arecP^. Here 
f (E) = [exp(E/k B T) + l}^ 1 is the Fermi distribution 
function at temperature T, and f = — r. The total charge 
current is obtained by summing over all contributions 
from the incoming waves, 



w = EE EE 7 - 



(28) 



=± "=± J=l,2 



Notice that / is also a function of A, e ex , T, Zj a and the 
geometry of the system (characterized by L and Wi/ 2 /s 
here) . Equipped with this form for the current across the 
interface, we are now in a position to explore it under 
various conditions. 



IV. RESULTS 

In this section we present major differences in trans- 
port properties resulting from SDIB, compared to the 
previous BTK studies^. The key physical quantity show- 
ing such differences is the normalized relative resistance, 
SR, defined as resistance difference between two cases of 
the device that have parallel (P) or anti-parallel (AP) 
magnetization alignment of the two F leads, 



SR 



Rap — Rp 



Iai 



(R AP +R P )/2 (7p+/ap)/2 : 



(29) 
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FIG. 2: (Color online) Relative resistance SR (defined as 
(Rap — Rp)/[(Rap + Rp)/2]) vs interface parameters Z. The 
solid curve is for a purely spin-dependent barrier (Zo = 0, 
Z s = Z) while the dashed curve is for a spin-independent one 
{Z = Z,Z S = 0). 




FIG. 3: (Color online) Contour plot of SR in the plane of 
Zq and Z 3 . The positive (lighter green) and negative (darker 
blue) regions are separated by a dashed contour indicating 
SR = 0. The solid contours are for SR = 1(T 4 , 10~ 3 and 



10 in the positive region as well as —10 
and — 10 _1 in the negative region. 



-10" 



-10" 



where /p(AP) is the total current obtained from Eq. (|28|) 
for the corresponding cases and -Rp(AP) = ^V-fp(AP) 



"he 



P (AP) case is characterized by the same (opposite) signs 
of the exchange energies in the two F leads, which leads 
to their majority spin species being the same (opposite). 
Previous calculations in Ref. 36 shows that SR is always 
non-negative, no matter how the barrier strength, ex- 
change energy, geometry of the device and temperature 
vary. This result is as expected of singlet-paired super- 
conductors and in and of itself would indicate that the AP 
case always carries more current than the P case given 
the same bias voltage, in conflict with the experimental 
results in Ref. 35, Here we show the exact manner in 
which the SDIB would alter this expected trend. 

In the following calculations, we take the majority spin 
in Fl to be t, ef = ±ef for the P (AP) case, and eV = 



() 



pi 




FIG. 4: (Color online) Relative resistance, 5R, vs exchange 
energy, e cx , at the interface parameters {Zo,Z„} = {0,0.7} 
(solid curve), {0,0.4} (dashed), and {0.7,0} (dotted). Inset: 
a zoom-out that shows 5R dropping to —2 at large imbalance. 



10 _2 A(at 0T) = 2 x lO~ 5 e F , which are typical values in 
experiments. We first consider zero temperature where 
the system is in the subgap regime and no quasi-particle 
transmission occurs. We analyze the effect of two factors, 
ferromagnetism as well as geometry, on which SR shows 
strong dependence and discuss how SDIB are relevant 
to the results. Then we study the finite temperature 
case where A varies in temperature and quasi-particle 
transmission contributes to the conductivity. 



A. Effects of ferromagnetism 

The effects of ferromagnetism on the charge transport 
emerge in two ways. The first is that the presence of e cx in 
the F leads results in a density imbalance in spin species 
and thus alters both the number of channels and the 
momentum of particles [see Eq. ( fl2| )] responsible for car- 
rying charge current. The second is that e cx in the inter- 
face region constitutes SDIB parameters Z (as discussed 
in Sec. [n]), which lead to different scattering phases be- 
tween majority and minority spins, and in turn to inter- 
ference in the scattering wave function, as we will see in 
subsequent discussions and results. 



Mathematically, this can be seen by noting in Eq. ( 24 ) 



that the current is a product of particle momentum, p, as- 
sociated with the first effect, and scattering magnitudes, 
\a\ and |6|, altered mainly by the second effe ct. The num- 

also reflects 



III) 



ber of channels, , fe (discussed in Sec. 
the imbalance effect and plays a role when we add up 
the probability current in Eq. (l27|) and (28). Both ef- 



fects contribute to a difference in behavior between the 
P and AP cases solely due to the presence of coherent 
crossed transport (CNR and CAR). This can be seen 
by noting that in the absence of these processes, for in- 
stance, for large separation between F leads, each lead 
can be treated independently and thus the P and AP 
cases would show the same results. On the other hand, 
the CNR and CAR processes distinguish the fact that in 
the P case the majority spin species are the same in both 
leads but are different in the AP case. We expect that 



the first effect is obscured in the low imbalance regime 
where the majority and minority species are less distin- 
guishable, while the second effect is suppressed when the 
interface barrier is nearly spin-independent (Z s ~ 0). 

In this subsection we choose the set of parameters 
Wi = W 2 = L = lOkp 1 , W s = 100fc F \ and zero tem- 
perature (T = 0). Figure [2] shows SR as a function of 
a purely SDIB (Z — and Z s varied, solid curve) and 
a spin-independent case [Z s = and Z varied, dashed 
curve) at a small imbalance of e GX = 0.0 lei?- We see 
that a purely SDIB causes an obvious variance of SR 
in both positive and negative values at an intermediate 
Z s . At large Z s the barrier is high enough such that 
most incidence is directly reflected while at small Z s the 
phase shift between majority and minority spins is close 
to zero, both of which make the difference between P 
and AP cases negligible and hence give SR ~ 0. For 
the spin-independent barrier, there is almost no differ- 
ence between P and AP, so SR is always flat and close to 
zero (still negative, reflecting the small imbalance effect 
as discussed before). A complete dependence of SR on 
both Z and Z s is shown in the contour plot in Fig. [3] 
We see that SR is negative for most of the parameter 
space and only becomes positive when Z s dominates. 

Although the maximum magnitude in the positive re- 
gion is one order smaller than that in the negative one, 
there is a wide enough parameter regime for 0.1% < R < 
1.5%, the same order magnitude observed in the exper- 
iment of Ref. 123 The most negative value is centered 
around the region of line Z = Z s where either majority 
or minority spins are subject to an almost transparent 
interface. This is because in the AP case different spin 
species in the two F leads see the transparent barrier, 
so the CAR is greatly enhanced and hence reduces the 
resistance compared to the P case. Notice that there is a 
small positive region at Z s < Zq, but the value is close to 
zero (SR < 10 -4 ) and becomes negative at Z s = 0. The 
contour plot remains unchanged if Z s — > — Z s , which 
reflects that only the magnitude of the spin-dependent 
component matters in the low imbalance regime. 

We now analyze the effect of density imbalance. Fig- 
ure [4] shows comparison of SR as a function of e cx for 
two purely SDIB conditions where SR is initially positive 
(Z s = 0.7, solid curve) and negative (Z s = 0.4, dashed) 
at zero exchange, as well as for a spin-independent condi- 
tion (Zq = 0.7, dotted). First we see that at large imbal- 
ance all curves drop to the highly negative region. This is 
due to a combination of number of channels and the mo- 
mentum for the minor species monotonically decreasing 
with the increase in e GX , as shown in Eq. (12). This de- 



crease reduces current carried by all Andreev processes 
involving one majority and one minority spin; CAR in 
AP case is not reduced as it involves two majority spins 
in the two F leads. At large imbalance in the P case most 
incoming current directly reflects and leaves small net 
current compared to AP case. Therefore SR drops and 
finally reaches its maximum value of —2 [see Eq. (29)] at 
which Ip = (see Inset). 
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FIG. 5: (Color Online) Local normal, local Andreev, crossed normal, cross Andreev currents, and their combined effect (blue 
long-dashed, green dotted, red short-dashed, purple dot-dashed, and brown solid curves, respectively) vs exchange energy, e cx . 
From top to bottom rows presents the current of the P case, 7p, AP case, Jap, and the difference, 51 = Ip — Tap, respectively. 
The first and second columns correspond to a spin-independent barrier and a purely SDIB of strength 0.7, respectively, while 
the third and fourth correspond to strength 0.4. The coincidence of two curves in some regions indicates the two local (crossed) 
currents almost equal to each other. 



The imbalance exercises a monotonic influence on the 
transport, explaining the curve for the case of spin- 
independent barrier. For an SDIB, however, SR ex- 
hibits an overall decreasing trend but a locally (e.g. in 
0.5 < e cx /ei? < 0.7) non-monotonic form as a function of 
e cx . We attribute this behavior to an interference effect 
of scattering through the SDIB. Differing interference ef- 
fects between the P and AP cases can thus either enhance 
or compete with the imbalance effect and dominate over 
it to make SR largely negative or even positive at low 
imbalance. 

We further investigate the imbalance effect by looking 
at four different contributions to the total current: the 
currents carried by (1) the incoming electrons and LNR, 
(2) CNR, (3) LAR, and (4) CAR [7 LN , I CN , I LA and 7 CA , 
respectively, defined in Eq. ( |B5| ) in Appendix [B] . Figure 
[5] shows these components and their combined effect as a 
function of e ox for P and AP cases as well as the difference 
between P and AP at various interface conditions. We 
first see that the SDIB cases have more current (lower 
resistance) than the spin-independent cases of the same 
barrier strength. For all cases the two local currents are 
the dominant contributions to the total current and are 
comparable in magnitude. Their monotonically decreas- 
ing trends are consistent with an increase in imbalance, 
except 7 LN in the low imbalance regime (e cx < O.llep) in 
the case of a purely SDIB of Z s = 0.4 (the left column), 
which reflects a stronger interference effect than the im- 
balance effect. The two crossed currents have such small 



contributions compared to the local current that they 
are almost indiscernible from zero in the figure, except 
I CA in the high imbalance regime in the two SDIB cases, 
which indicates a great enhancement in long range en- 
tanglement. The difference in current between P and AP 
Ip — Iap [directly related to SR via Eq. ( 29 )] , 



SI 



is one to two orders smaller than either Ip or Jap arid is 
highly sensitive to the value of the exchange field. 

The two spin-independent cases have a similar behav- 
ior: at low imbalance the compensation of dominant con- 
tributions of (5/ LN and SI makes the total SI close to 
zero (still negative), while at high imbalance i5/ CA has 
a large negative contribution and SI hA is suppressed, as 
supported by the imbalance effect. 

The purely SDIB cases have various dominant contri- 
butions at low imbalance. At Z s — 0.7, the dominant 
contribution is positive <5/ LN , making the sum of the con- 
tributions positive, while at Z s = 0.4 the strong negative 
SI CN makes the sum negative. These results can only 
be attributed to the interference of scattering through an 
SDIB. At high imbalance, SI CA dominates as in the spin- 
independent case but is much larger, due to the greater 
enhancement of I CA in the AP case. Notice that all the 
curves have kinks at the same positions - this is due to 
the reduction in the number of scattering channels by one 
at each kink. 

In brief, we analyzed the effects of spin imbalance and 
scattering through an SDIB, both induced by ferromag- 
netism. The former monotonically lowers SR toward neg- 
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FIG. 6: (Color online) Relative resistance, 8R, vs the separa- 
tion between two F leads, L, with fixed Wim (Wi = Wi = 
lOfc^ 1 ) at the interface parameters {Zq,Z s } = {0,0.7} (solid 
curve), {0,0.4} (dashed), and {0.7,0} (dotted), same as pre- 
sented in Fig. [4] Relative distance from a reference are marked 
on top of graph, showing the oscillatory period of n. 



ative values and dominates in the high imbalance regime, 
while the latter exhibits its influence in positive or nega- 
tive directions, where the trend is revealed in the scatter- 
ing magnitudes. A system with an SDIB can sustain a 
significant large CAR current in the AP case, reminiscent 
of a great enhancement of long range entanglement. 



B. Effects of geometry 

The signature of interference due to scattering be- 
tween the two F leads implies a strong dependence of 
the transport properties on their geometry. We exam- 
ine SR as a function of the widths, W\ and W%, and 
the separation, L, of the leads at a low imbalance value 
of e GX = O.Qlep where the scattering effects dominate. 



Here we fix Ws = lOOfc^ . Figure [6j shows SR vs L 
at three barrier conditions, {Z ,Z S } = {0,0.7}, {0,0.4}, 
and {0.7,0}, as presented in Fig. [4] For the two SDIB 
cases, the curves modulate and exponentially decay to 
zero at large L. The oscillation period is 7r (in units of 
kp 1 ), which results from a 2tt modulation in the scat- 
tering amplitudes, which in turn owe to the oscillatory 
behavior of matrix O in Eq. (22). Since the quantity A 

and L, sum- 



in Eq. ( 19 ) is a sinusoidal function in W± 



/2/S 



ming the product of A in calculating f2 directly indicates 
a relation between the interference and the geometry of 
the system. The fast decay is due to the lack of corre- 
lations between the F leads at separations larger than 
the superconducting coherence length, where scattering 
hardly shows a difference between A and AP cases. For 
the spin-independent barrier, the curve is flat and close 
to zero as expected. Notice that the modulations in the 
current for P and AP cases still exists (as discussed in 
Ref. l36|) . However, they modulate in-phase and cancel 
out in the form of SR. 

Figure [7] shows contour plots of SR in the plane of 
W 1 and W 2 at L = (Wi +W 2 )/2 for the pure SDIB 
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FIG. 7: (Color online) Contour plot of 5R in the plane of 
Wi and Wi. Data presented for a pure SDIB of Z s — 0.7 
(Panel a) and 0.4 (Panel b). The lowest contour is zero, be- 
low which the value is negative (corresponding to the black 
region). The magnitude difference between the contours is 
0.2%. Relative distance from a reference are marked on top 
and right of graph, showing the oscillatory period of n. 



cases of Z s = 0.7 (Panel a) and 0.4 (Panel b). We see 
SR modulates with a period of ir in both W\ and W2 
directions, for the same reason as for variation in L. In 
(a), SR is always positive and becomes close to zero at 
large W\ and W 2 . In (b), the value is mainly positive 
with a slower decay as W\ and W% go large. Both plots 
show a wide parameter range for positive SR in contrast 
to the small negative value expected in the case of spin- 
independent barriers. 

Finally we keep relative ratios between Wi, W 2 , Ws, 
and L unchanged and enlarge the whole device (making 
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FIG. 8: (Color online) Relative resistance, 8R, vs pure SDIB 
Z s . The solid curve (axis on the left of the graph) presents 
a purely SDIB system of Wi = W 2 = L = 10k' 1 and Ws = 
lOOfc^ 1 (same as presented in Fig. [2]). The dashed and dotted 
curves (axis on the right of the graph) are for devices with 
sizes 10 and 30 times that for the solid curve, respectively. 



FIG. 10: (Color online) Relative resistance, SR, vs the tem- 
perature, T, at the interface parameters {Zo,Z a } = {0,0.7} 
(solid curve), {0,0.4} (dashed), and {0.7,0} (dotted), con- 
vention as presented in Fig. [9] Notice that the dotted curve 
almost coincide with the axis. 
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FIG. 9: (Color online) Charge current, I, vs the tempera- 
ture, T, for the P case at the interface parameters {Zo, Z s } = 
{0,0.7} (solid curve), {0,0.4} (dashed), and {0.7,0} (dotted), 
as presented in Fig. [4] The star signs indicate position of the 
maximum point of the corresponding curves. 

it closer to some of the realistic systems^}). Figure [8] 
shows SR as a function of a pure SDIB parameter Z s 
for three difference scales of the device: the one of W\ = 
W 2 = L = lOfc^ 1 and Ws = lOOfc^ 1 (solid curve, same as 
presented in Fig. [2]) as well as 10 and 30 times the case of 
the solid curve (dashed and dotted curves, respectively). 
When the system becomes an order of magnitude larger, 
the positive SR region expands but its value is more than 
one order smaller. Because larger sizes result in more 
conducting channels and hence more interference among 
them, this trend can smear out the difference between 
P and AP cases and leave a small SR. To improve the 
modeling for a real system that can still have SR of a 
few percent, new ingredients such as the dependence of 
Z parameters on spatial coordinate or on channels would 
be needed. 



C. Effects of temperature 

Finally, we briefly study finite temperature effects. 
Here, two factors need to be incorporated: (1) quasi- 



particle (hole) transmissions (QPT and QHT, respec- 
tively) that arise from thermal fluctuations and partic- 
ipate in charge transport as well as (2) the decrease in 
the superconducting gap as a function of temperature, 
A(T) = A(0)v/1 - T/T c . Here, we focus on the effect 
of the SDIB and thus choose the case of low polariza- 
tion for which the effects of SDIB are most prominent. 
We plot the current in the P case, Jap vs temperature, 
T, for different barrier conditions in Fig. [9] All cases 
show suppressed conductance for a significant tempera- 
ture regime close to T c due to the vanishing superconduc- 
tivity. The spin independent interface shows local peaks 
which we attribute to thermal population of new chan- 
nels. Spin dependence of the barrier suppresses these 
peaks. We also mention that in accordance with BTK 
expectations^, while the SDIB curves seem to decrease 
as a function of temperatures, the maximum (star signs 
in graph) occurs at non-zero but very low temperature, 
indicating an re-entrant effect. 

To characterize the difference between P and AP cases, 
we plot SR vs T for the same set of parameters in Fig. [10] 
We see that SR can be smooth, widely fluctuating, or 
even endowed with a sign change for different SDIB con- 
ditions. The salient feature here is that in the spin inde- 
pedent case, due to the low polarization, SR is exceed- 
ingly small and constant on the scale shown. However 
the SDIB cases show that temperature can alter the sign 
of SR and what might have started as negative at low 
temperature can reverse the trend. The lack of trend has 
also been found in different samples in experiment!^. 

V. CONCLUSION 

In this paper, we studied electronic transport physics 
in superconductor-double ferromagnet (S-FF) junctions 
with spin-dependent interface barriers (SDIB) using the 
Blonder- Tinkham-Klapwijk (BTK) treatment. We pro- 
posed spin-dependent Z parameters modeling SDIB for 
the BTK calculations and discussed its relation to micro- 
scopic physical quantities at the interface. Our extended 
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model is one of the first incorporating SDIB for S-FF 
junctions and shows that SDIB can increase crossed An- 
dreev reflection (CAR) current. We found that SDIB can 
cause interference in wave functions of CAR and other 
scattering processes that can either collaborate or com- 
pete with the imbalance effect in the F leads. This com- 
petition can cause more current carried by F leads of 
parallel (P) magnetization than by the antiparallel (AP) 
case (resulting in a positive SR), which does not occur in 
the case of spin-independent barriers. Our calculations 
showed that positive SR is stable in a sufficiently large 
parameter regime as a function of Z parameters, mag- 
netization, geometry of the system and temperature. At 
low magnetization, the effects of SDIB are responsible for 
the different transport between P and AP cases, result- 
ing in positive or negative SR that is sensitive to the Z 
parameters. At high magnetization, 5R is always nega- 
tive due to the population imbalance in the F leads, and 
SDIB contributes large additional enhancement of CAR 
in the AP case. Variance in the geometry of the junction 
results in oscillations in 5R, which reveal the interference 
effects on the transport physics. 

Our studies suitably describe recent experiments of 
Ref. |35] and can be used as a guide for further investi- 
gation in such settings, including the possibility of con- 
trolled tuning of S-F interface properties in the presence 
of artificially grown SDIB barriers. Our model easily 
extends to a wide class of multi-terminal hybrid junc- 
tions. An important ingredient to take into account in 
the future, particularly to explain differing P versus AP 
behavior, is to incorporate induced spin-triplet correla- 
tions together with the spin-flip scatterer d 17 ! 42 * ^. An 
immediate extension of our studies would be to evaluate 
This would be particularly relevant to currently intense 
explorations of coherent properties and long range entan- 
glement in superconductor heterostructures. 



injected spin passes through the barrier region and accu- 
mulates a phase shift <p in the interface region, compared 
to the case without a barrier. In a semi-classical pic- 
ture, the phase shift due to the variance of the exchange 
energy, Se ex (x), is 



<(*(*)) 



dt = 



dx 



_ i{5e 



flVF 



(Al) 



where io is the total time to pass the barrier, vp is the 
Fermi velocity, Se e *(x(t)) is the local variance of exchange 
energy when the particle passes the position x at time t, 
and (<5e ex ) is its average over the width of the interface, 
£. Here x = vpt given that the energy variance is small 
compared to the Fermi energy. 

Now we compare the phase shifts with those of a scat- 
tering problem having the (^-function potential, 

-~d 2 x + hv F Z s 5(x). (A2) 
2m 

By solving the Schrodinger equation for an incoming 
wave around the Fermi surface, we obtain Z s as a func- 
tion of the phase shift of the transmitted wave, 

Z s = - tan 0. (A3) 

Substituting Eq. (A3 ) into Eq. ( Al ), we derive the aimed 
relation of Eq. ([8| . The tangent function implies that Z s 
in non-monotonic in and sensitive to the exchange energy 
and microscopic length. Notice that the majority and 
minority spins have opposite exchange energy which leads 
to different signs for the spin-dependent Z parameter. In 
this paper we choose zhZ s for majority (minority) spins, 
as shown in Eq. Q. 
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Appendix B: I-V relation 

In this section we detail the calculations for current- 
voltage relations in t he SF F junction following the stan- 
dard BTK formalisrrpJ^SI anc j identify the current con- 
tributions for each scattering process. The total current 
contributed by an incoming wave characterized by indices 
{rajn} is given by 



I TO 



At / dE[U(E)-f^(E)], 



(Bl) 



Appendix A: Semiclassical picture for SDIB 

In this section we use a semiclassical picture mod- 
eling the spin-dependent component of the SDIB pa- 
rameter, Z s , in Eq. |7|). We derive Z s by relating the 
phase shifts between scattering through the phenomeno- 
logical (5-function barrier and that through the micro- 
scopic Zeeman potential at the interface, provided the 
spin-independent component Z — 0. Microscopically, an 



where A is a constant associated with density of states, 
Fermi velocity as well as an effective cross-sectional area 
and /->(<-) denotes the incoming (outgoing) populations 

U(E) = fo(E-reV), 

fUE) = 



E 

->'}=±{t 

EE! 

r'—zt n' 



EE 



Jt a j n 
rcrjn 



,<t} 3' 



(B2) 

fo(E-r'eV) 



Jr'Jjn fo(E). 



(B3) 
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Here we replace the distribution function for holes in the 
F regions as l-fo(-E-eV) = f (E+eV). By substitut- 
ing Eqs. (|B2| and (IBS! into Eq. (|B1|) as well as applying 



a conservation relation for the probability currents, 



E EE 



Tajn 



EE 



jSyr'n' 
rcrnj 



{t',ct'}=±{t,<t} j' 



= 1 



we eliminate the contributions from quasi-particle trans- 



missions and hence derive the I-V relation of Eq. (27 1, 



as a function of only processes in F regions. By summing 
over the incoming degrees of freedom, we obtain the total 
current as in Eq. (28), or simplified at zero temperature 
as 



<t=± n i=l,2 ' 



n' j'=l 



n j 
nj 



J 



-\-<rn / j' 



(B4) 



From Eq. (B4) we identify four components of the to- 



tal current in the subgap regime below: the local nor- 
mal current, 7 LN , contributed by the incoming wave and 
LNR, as well as the crossed normal current, J CN , the local 
Andreev current, J LA , and the crossed Andreev current, 
7 CA , contributed by CNR, LAR, and CAR, respectively. 



rLN 



jCN _ 



jLA(CA) 



eV 



J 



A EEE 

cr=± n,n' j—1,2 q 
a— ± n,n' j— 1,2 q 

•'EE E J^\rz/ G) 



+crn j 
anj 



J +crnj 
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